! 1D - function integration
module func_int_mod; 
use test_mod;
use grid_mod;
use vec_mod;
implicit none;
logical, save :: FUNC_INT_TRACE_ENTRY = .TRUE.;
logical, save :: CACL_DRV_P5_TRACE = .TRUE.;
real(R14), dimension(4), save :: INT_P5 = & 
    (/14D0, 32D0, 12D0,  32D0 /);

contains 

real(R14) function calc_int(xq, pf) result(res);    implicit none;  
  real(R14), dimension(:), pointer :: xq, pf;      
  call dbg_info("-->calc_int()", FUNC_INT_TRACE_ENTRY);
  res = vec_dot_prod(xq, pf);
  return;
end function

real(R14) function calc_int2(xq, pf, pf2) result(res);    implicit none;  
  real(R14), dimension(:), pointer :: xq, pf, pf2;      
  call dbg_info("-->calc_int2()", FUNC_INT_TRACE_ENTRY);
  res = vec_dot_prod3(xq, pf, pf2);
  return;
end function

real(R14) function calc_int3(xq, pf, pf2, pf3) result(res);    implicit none;  
  real(R14), dimension(:), pointer :: xq, pf, pf2, pf3;      
  call dbg_info("-->calc_int3()", FUNC_INT_TRACE_ENTRY);
  res = vec_dot_prod4(xq, pf, pf2, pf3);
  return;
end function

! new integration quadrature for px-step-grid
function new_step_qdr(px) result(res);    implicit none;  
  real(R14), dimension(:), pointer :: px, res;      
  call dbg_info("-->new_step_qdr()", FUNC_INT_TRACE_ENTRY);
  res => new_vec_from(px);
  call step_qdr_p5(px, res);
  return;
end function

! 5-point integration quadratures
! Not tested for dimension(J:J2), where J != 1
SUBROUTINE step_qdr_p5(px, res);    implicit none;  
  real(R14), dimension(:), pointer :: px, res;      
  real(R14) stx, h2; ! step x
  integer xn, mini, maxi, i, k;
  real(R14), dimension(4) :: wa; ! work array

  call dbg_info("-->step_qdr_p5()", FUNC_INT_TRACE_ENTRY);
  xn = SIZE(px, 1);  mini = LBOUND(px, 1); maxi = UBOUND(px, 1);
  call assert_true("xn >= 5", xn >= 5);
  call assert_true("must be n=k*4+1; MOD(xn-1,4)==0", &
    MOD(xn-1, 4) == 0); ! 
  call assert_true("xn == SIZE(res, 1)", xn == SIZE(res, 1));
  call assert_true("mini == LBOUND(res, 1)", mini == LBOUND(res, 1));
  call assert_step_grid_4(px);
  stx = px(mini+1) - px(mini); ! step
  h2 = D2 * stx / 45D0;
  wa = h2 * INT_p5; ! load weights for the given x-grid
 
  call dbg_vec("px=", px);
  ! not tested but trying to have independent mini
  do i = mini, maxi
    res(i) = wa(1 + MOD(i - mini, 4)); 
  end do;
  res(mini) = res(mini) * D05;
  res(maxi) = res(maxi) * D05;  
  return;
end SUBROUTINE

end module 